Simulating Van der Waals- interact ions in 

"S" 

O water/hydrocarbon-based complex fluids 

(N 

> 

O 

00 

' (-('! I. Pasichnyk 

O 

a , 

^ . Max Planck Institute for the Physics of Complex Systems, Noethnitzer Str. 

^\ 

^— > ! 38, D-01187 Dresden, Germany 

R. Everaers 

^' 

^ ■ Universite de Lyon, Laboratoire de Physique, Ecole Normale Superieure de 

Lyon, CNRS UMR 5672, 46 allee d'ltalie, 69364 Lyon Cedex 07, France 



> 

OO . 

Max Planck Institute for the Physics of Complex Systems, Noethnitzer Str. 

' ■ 

38, D-01187 Dresden, Germany 

T— I ■ 

O ■ A.C. Maggs 

> 

U ■ Laboratoire de Physico-Chimie Theorique, UMR CNRS-ESPCI 7083, 10 

rue Vauquelin, F- 75231 Paris Cedex 05, France 



Abstract 

In systems composed of water and hydrocarbons Van der Waals- 
interactions are dominated by the non-retarded, classical (Keesom) 
part of the Lifshitz-interaction; the interaction is screened by salt and 
extends over mesoscopic distances of the order of the size of the (mi- 
cellar) constituents of complex fluids. We show that these interactions 
are included intrinsically in a recently introduced local Monte Carlo 
algorithm for simulating electrostatic interactions between charges in 
the presence of non-homogeneous dielectric media. 

1 Introduction 

In the special case of water-hydrocarbon systems, which notably include 
biological systems, the weak optical contrast between water and many hy- 
drocarbons leads to Van der Waals interactions which are dominated by the 
classical (Keesom)-contribution [T]. Within the Lifshitz formalism it is possi- 
ble to perform analytical calculations only for continuum descriptions of sim- 
ple geometries such as planar interfaces and lamellar structures [21 [3l HI [5] . In 
the opposite extreme of atomistic Molecular Dynamics simulations, the rele- 
vant partial charges on the water (solvent) molecules are treated explicitely. 
This results in proper treatment of the Keesom contribution of the Van der 
Waals interaction (because in the microscopic simulation the microscopic 
dipoles are fluctuating). However, current all-atom simulations are limited 
to the nanosecond timescale, while the physical processes can take much 
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longer. Coarse-grained descriptions using implicit solvent models can help 
to close this gap [6]. Typically, electrostatic interactions are calculated from 
a (macroscopic) dielectric theory [7] while Van der Waals forces are mod- 
elled via (effective) Lennard- Jones interaction with a cutoff of the order of 
the "grain size" a of the coarse-grained system [8]. However, neglecting 
the collective origin of the Van der Waals interactions, namely electrostatic 
interactions between fluctuating charge distributions, gives rise to several 
problems: 

- in complex fluids Van-der- Waals interactions extend over distances 
which are comparable to the size of interacting objects [9], [lO]. For 
polymers of amphiphilic systems organizing into micellar or lamellar 
structures, the characteristic length scale can be much larger than the 
size of constituting monomers. 

- the results obtained by the approximated implicit solvent models are 
very sensitive to the value used for the dielectric constant, which turns 
out not to be a universal constant but simply a parameter that depends 
on the model used ITTI: 



the effect of screening by salt of a classical contribution to the Van-der- 
Waals interaction [12] is not included. 



Recent publications [131 [HI [IS] reformulated the problem of electrostatic 
interactions between charges in the presence of non-homogeneous dielectric 



media: Long ranged electrostatic interactions between charges are gener- 
ated dynamically via local interactions between charges, the medium and 
the electric field. Previously in [H] the interaction between small dielectric 
inhomogeneities was considered in the dilute limit. Here we generalize the 
proof to arbitrary densitites and show that the method implicitly generates 
the many-body effects in the zero frequency part of the Lifshitz interaction 
regardless of the system density. 

The paper is structured as follows: After a short review of Lifshitz theory, 
in Section [3] we present the central result of the paper - the theoretical basis 
of the simulation method and its relation to Lifshitz theory. To be able to 
treat systems with general geometries we have to validate our method for 
the case of simple systems where the theoretical result can be used for the 
comparison. Therefore we have chosen a triple slab system since for this par- 
ticular geometry one can develop and test a technique for correct simulation 
and thermodynamic integration. In Sec. H] we present our simulations and 
compare to analytic theory for the triple slab geometry. 

2 Theoretical background 

2.1 Lifshitz theory 



Dzyaloshinskii et al [12] recasted Van der Waals forces in terms of interac- 
tions between continuous dielectric media, mediated by the electromagnetic 
field. The result corresponds to summing a series of multi-body interactions 



between fluctuating charges. In describing Van der Waals interactions the 
specificity of the condensed medium is completely taken into account by us- 
ing its dielectric function e(u;), u being the frequency of electromagnetic field. 
In particular, the electromagnetic field fiuctuation free energy JF is given by 



T=kBTY^\nV{iin) 



(1) 



n=0 



where ks is the Boltzmann constant and T is the temperature. The n sum- 
mation is over bosonic Matsubara frequencies ^„ = 27rnkBT/h. The prime in 
the summation refiects the fact that the n = term is taken with a weight 
1/2. The secular mode equation (or dispersion equation) V{i^n) = gives the 
eigenfrequencies of the electrodynamic field modes in the specified geometry. 
For the case of two plane parallel half-spaces with dielectric constant ei 
separated by the gap of length / and dielectric constant €2 one can explicitely 
derive the free energy (per unit area) of the interaction [12] : 
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The susceptibilities e = ^{i^n) are evaluated on the imaginary frequency axes. 

The fluctuation driven electromagnetic interactions may be classical or 
quantum in origin. Usually the temperatures of interest for condensed media 
are low compared with hujQ, where ujq is a typical frequency in the system 
which is often in the ultraviolet (Tq ~ huJo/kB ~ 7 x lO^K) [17]. In most 
condensed matter systems Van der Waals interactions are thus dominated 
by quantum fluctuations. Important exceptions occur in mixtures of polar 
liquids (e.g. water) and hydrocarbon based (macro)molecules, a situation 
of considerable interest for biological and biophysical problems. This is a 
consequence of two effects. Firstly, there is low contrast between the dielectric 
response of water and hydrocarbons in the optical part of the spectrum. 
Secondly, there is a large contrast at /oty-frequencies due to orientational 
fluctuations of dipoles in polar liquids. 

If one works in the gas phase, rather than in condensed media, and con- 
siders the interaction energy between two water molecules in vacuum, the 
corresponding classical Keesom forces [9j at room temperature are charac- 
terized by the prefactor to the interaction in 1/r^: C^'^^'^°'^ = 96 x 10~^^Jm^ 
considerably larger than the quantum (known as dispersion) contribution 
with the strength C^^^^ = 33 x 10""^^ J m^ [H]. As a result the zero-frequency 
contribution to the Van der Waals interaction in water-hydrocarbon systems 
dominates and gives approximately 60% of the net interaction potential [9]. 

When one drops, in the sum Eq. [21 the terms for which n 7^ (which are 



essentially quantum mechanical) one finds 



oo , 2 



"^ 

Eq. [5]can be derived using a different approach OH]. Dean et al. [1] have 
shown that if one considers a thermal field theory for the field i\} with purely 
electrostatic Lagrangian 



l^\iA = \ I eir)iV^fdh (6) 



the zero frequency Lifshitz term can be obtained from the partition function 
of field V': 

Z= [ d[^|J]exp{(3C[^P]) (7) 



where (3 = l/{kBT). After changing in the latter formula the axes of func- 
tional integration via ip — > —i(f) one recovers the partition function of the 
dielectric system |1]: 

Z= f d[(p]exp(- f (P\/e{v)\/(pdhj 
= [det(-Ve(r)V)]-^/' 

where det(— Ve(r)V) is formally understood as the product of eigenvalues of 
operator — Ve(r)V [19]. Finally Eq. [5]can be recovered for the case of planar 
geometry if one calculates the free energy from the usual thermodynamic 
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relation T = —ksTlnZ. 

In all that follows we will consider only the zero frequency term (see Eq. [5]) 
of Lifshitz interaction and will drop the subscript ra = 0. Furthermore we 
will consider only static susceptibilities and will not specify the frequency 
argument of e. 

2.2 Triple slab geometry 

A triple-slab geometry (Fig. [1]) belongs also to the class of analytically solv- 
able geometries. We will use it to compare our simulations to known results. 
It is also easy to treat in periodic boundary conditions. As the free energy 
is an extensive quantity [20], it contains a volume contribution as well as 
a surface contribution. We are interested in the distance dependence of the 
surface part of free energy. The triple-slab geometry allows one to change the 
distance between two slabs without any changes in the volume of dielectric 
materials in finite systems. Hence the volume contribution in such a system 
can be easily separated from the surface free energy we are interested in (see 
SecS]). 

One considers two slabs of material with dielectric constant ei of thickness 
h and area L^ which are separated by a distance h. The dielectric constant 
of the external medium is given by €2 (see Fig. [1]), so that 



e{z) = e2 + (ei - e2)d{z) + (es - ei)e{z - b) 

(9) 
+ (ei - e2)9{z -h-h) + {e2- e,)9{z - 2b - h) 
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Figure 1: Symmetric triple layer. 



where d{^z) is the Heaviside step function. 

The interaction energy per unit area is given by [3l H] : 



where 



T{h,h) 



Air 



dp pin 



A^fl 



g-2fep\2g-2hp 



(1 - A2e-2bp)2 



£2 -ei 



(10) 



(11) 



62 + ei 

One can also consider the pairwise approximation to the general result 
given by Eq. [TOl The major contribution to the integral with respect to p 
comes from the saddle point p ^ 0. Hence we can write, since A^ < 1, 



ln(l - A^e-'^P) ~ -A^-^'^P, 



(12) 



and carry out the integration with respect to p to find [3] 






In the case of two infinite slabs 6 ^ oo we have the usual result of Hamaker 
theory [21]: 

where A is the classical part of the Hamaker constant. 



3 A local Monte Carlo algorithm for gener- 
ating Van der Waals interactions 

In the following we describe our simulation method. At zero temperature 
the Coulomb interaction results from minimizing the energy 



1 /■ D^ 
^ = -2 7^"'' (15) 



where e(r) is assumed isotropic and D is the electric displacement constrained 
by Gauss law, V ■ D = p; p is the external charge density. We assume 
throughout the paper that the dielectric constant of the vacuum is eo = 1. 
This constrained minimization problem for D can be solved with the help of 
a Lagrange multiplier 0(r) by looking for stationary points of the functional 



D _, 

W[D] = y"|^-0(r)(V-D(r)-p(r))|d=^r (16) 

and is given by 

U, = \je{v){Vcl^fdh (17) 

is the solution of the Poisson equation 

V ■ (e(r)V</)) = -p (18) 

The true interest of the formulation appears in Monte Carlo since one 
can generate local dynamic systems which sample the partition function 

r ^ 
Z= Yl d\ Y[ I?D(r)(5(V ■ D - p(r))e- ^ /°V<r)d3r ^^g^ 

•^ i=l r 

Following [13] we discretize the system placing particles on a simple cubic 
lattice with vector fields such as D on the links. This formulation is numer- 
ically efficient because a local variation in p requires only a local update of 
the field D. For problems involving dielectric inhomogeneities (macropar- 
ticles with dielectric constant differing from the surroundings) one has to 
choose an appropriately interpolated value of the dielectric function. The 
dielectric function is placed also on the link (Ref. [H]) and is given by the 
harmonic average 

A.l,^ (20) 



^nfi ^n *-Ji+/i 
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where {n/i} is the hnk which goes from the site n in the //-direction, // 
1, 2, 3. In the absence of charged particles the partition function Eq. 
becomes 

Z= /"jJPD(r),5(V-D)e-2/D'AW'^''- (21) 

•^ r 

Introducing an auxiliary field to implement the Gauss' law constraint and 
using the identity 2ti5{x) = J e"^"^ d(p the last equation is equivalent to 

•^ r r 



c^u<^r m 



(22) 



3/2 / TT^^.-275/^W(V</')^'i^ 



e 2/3. 



= C3 [det(-V • e(r)V)]-^/' 

where the constants Ci, C2 and C3 are of no further interest [1|. Comparing 
Eq. [22] and Eq. [8] we conclude that our method produces the zero-frequency 
term of the Lifshitz interaction. One has to note, in spite of the fact that 
intermediate expressions in deriving Eq. [22] contain complex contributions, 
the algorithm directly samples the real constrained partition function given 
by Eq. M 



^In the case of moving particles the variations of the term ^/e can add nontrivial 
contributions to the contact energy. 
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4 Numerical validation 

We simulate a triple slab system Fig. [T] with two values of dielectric constant 
of the media: ei = 2 and ei = 78 using different values of the lattice constant 
a. The dielectric constant of the intermediate region is e2 = 1 for both cases. 
The size of the box is L = 15.0, the width of the slab h = 1.0. 

In order to calculate the free energy we perform a thermodynamic inte- 
gration [23]. For the reference system (denoted by /) the uniform dielectric 
constant €2 = 1.0 is assigned. We sample the system with the potential energy 
U which depends linearly on the coupling parameter A: 



W(A) = (1 - Xpi + XUji 



[1-X) I d'r— + X I dh^ 
' ' 262 J 2e(r) 



(23) 



For A = 1 we recover our system of interest (denoted by //). The system 

with energy U{X) is equivalent to the system with the following dielectric 

function: 

6.(A,r)= ,,J['^'' ,,, (24) 

e(r) + A(e2-e(r)) 

Finally, the free energy difference between systems II and I can be found 
from the following expression: 



r-(i) 



where (...);!^ denotes an ensemble average for a system with energy Eq. 
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The numerical calculation of a free energy is always demanding. We have 
approximated the integral in Eq. [25] by a summation over 20 intervals in A. 
The fluctuating field D is sampled by a worm algorithm [211 125] • Simulation 
at each A point involved an equilibration period of 800 sweeps, where a 
sweep is 20 worms, followed by a consequent run of another 800 sweeps 
configurations. The error bars and average values of free energy have been 
calculated from 500 values of free energy. Simulations were performed on an 
AMD Opteron 2.4GHz processor. Total simulation time for a one measured 
point was 2 days for a = 1.0 and 24 days for a = 0.5. 

The free energy calculated in this way gives the full contribution which 
includes self-energies of individual slabs as well as the interaction energy 
between slabs. In contrast Eqs. [T3l and [TO] represent only the interaction part 
of the excess free energy. In a system with periodic boundary conditions it is 
difficult to calculate the limit h ^ oo which corresponds to calculating the 
self-energy part. Therefore we perform an interpolation of our simulation 
results by the function Eq. [13] and extrapolate them to the region h ^ oo 
to find the asymptotic value. Further we subtract this contribution from the 
total free energy Eq. [25] Of course such a procedure leads to small deviations 
from the analytic curve which can be clearly seen on the corresponding plots. 
We are interested in observing the free energy of the system as we change the 
separation between slabs. Our results are shown in Fig. [2] for ei = 2.0 and in 
Fig. [3] for ei = 78.0. In both cases, a comparison of results for two different 
values of the lattice constant a = 0.5; 1.0 shows that the errors due to lattice 
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Figure 2: Free energy of the slab at ei = 2. 



10 



d 10° 



Xi 

j^ 



10" 



10' 



1 ' ' 


■ 


•N^^^ 


A a=l 


" 


- ^^^^^ 


D a=0.5 


^ 


^^r^\^A 


— analytic result 


: 


^v^^ 


— pairwise 


- 


N^ ^**Ss»^^ 








"S, ^*^^ 




•Sn ffl***Sfc. 








Xj ^»jR 


: 


V 1 p** 


- 


p 




N 


1. T T 




•s . 


f^ 


L 


[ 


] -r 


■ 




•s ■ 


'Nr 


1 








V 


B 


J^ 






j- 




s 


N 




—_ 


: 




N 


J 


k 


■ ■ j 


: 


. 


. 




hi 


A i.- 


■ 




N 


\ 


■ 




■» 


sh 


1 1 1 1 1 1 1 1 


1 


1 


1 



Figure 3: Free energy of the slab at ei = 78. 



14 



discretization are small. In particular, the data reproduce the analytic result 
Eq. [To] which differs significantly from the pairwise (Hamaker) curve Eq. [13] 
for large dielectric contrasts. 

5 Conclusion 

We have shown that a recent Monte Carlo algorithm for the simulation of elec- 
trostatic interactions in heterogeneous dielectric media implicitly generates 
the zero-frequency part of the Lifshitz interaction including all many-body 
effects. The interactions make the dominant contribution to the Van der 
Waals attraction in hydrocarbon-water systems as they are typically found 
in soft and biological condensed matter [9]. The method is easily applica- 
ble to systems with interfaces and spatially varying dielectric constants of 
arbitrary geometry and allows the inclusion of fixed and free charges. 
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